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A cascade model is described based on multiplier distributions determined from 3D direct numer- 
ical simulations (DNS) of turbulent particle laden flows, which include two-way coupling between 
the phases at global mass loadings equal to unity. The governing Eulerian equations are solved 
using psuedo-spectral methods on up to 512 3 computional grid points. DNS results for particle 
concentration and enstrophy at Taylor microscale Reynolds numbers in the range 34 - 170 were 
used to directly determine multiplier distributions on spatial scales 3 times the Kolmogorov length 
scale. The multiplier probability distribution functions (PDFs) are well characterized by the (3 dis- 
tribution function. The width of the PDFs, which is a measure of intermittency, decreases with 
increasing mass loading within the local region where the multipliers are measured. The functional 
form of this dependence is not sensitive to Reynolds numbers in the range considered. A partition 
correlation probability is included in the cascade model to account for the observed spatial anti- 
correlation between particle concentration and enstrophy. Joint probability distribution functions 
of concentration and enstrophy generated using the cascade model are shown to be in excellent 
agreement with those derived directly from our 3D simulations. Probabilities predicted by the cas- 
cade model are presented at Reynolds numbers well beyond what is achievable by direct simulation. 
These results clearly indicate that particle mass loading significantly reduces the probabilities of 
high particle concentration and enstrophy relative to those resulting from unloaded runs. Particle 
mass density appears to reach a limit at around 100 times the gas density. This approach has 
promise for significant computational savings in certain applications. 

PACS numbers: 47.61. Jd, 47.27.E-, 47.27.eb 
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I. INTRODUCTION 

The study of turbulent flows incorporating heavy par- 
ticles in suspension (particles with finite stopping times) 
is an important endeavor that has both fundamental and 
practical relevance to many scientific and engineering 
problems. Such flows have been investigated mainly in 
numerical simulations where detailed statistical analysis 
of the flow fields is possible [H, 0> S B| These simulations, 
limited to relatively low Taylor microscale Reynolds num- 
bers Re\ (~ 40), demonstrated that particles whose fluid 
response times are comparable to the lifetime of the 
smallest turbulent eddies produce a highly nonuniform 
field with intense regions of concentration. Preliminary 
indications were that the feedback from such concentra- 
tions of particles could locally damp turbulence - how- 
ever, the role of this "mass loading" effect in determining 
the statistical distributions of particle density and vari- 
ous fluid scalars has not been thoroughly studied. Ex- 
perimental investigations of turbulence modification by 
particles have demonstrated that the degree of turbu- 
lence damping increases with particle mass loading and 
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concentration [4J. 

The phenomenon known as intermittency can be de- 
scribed as intense fluctuations, on small spatial and tem- 
poral scales in the turbulent field, that contribute to 
the exponential tails of probability distribution functions 
(PDFs) of scalars such as velocity increments and gradi- 
ents J3S,0|> dissipation Q, pressure [H, [l(|, enstrophy 
[ill fl21 | and velocity circulation [13j |. Intermittency in 
the density field of preferentially concentrated particles 
has also been observed and studied [3, EH • 

Although intermittency in turbulence still lacks a com- 
plete theoretical understanding, progress has been made 
with phenomenological models that capture intermit- 
tency in a cascade process. Richardson [16[ and later Kol- 
mogorov [l7| suggested that such models might be used 
to explain the process of eddy fragmentation initiated 
by unstable large scale structures in a turbulent fluid. 
Intermittency in the context of fragmentation though a 
cascading process has been studied for large-scale gravi- 
tating masses [l8| and velocity increments in turbulence 
Simple cascade models were explored by Meneveau 
and Sreenivasan l20ll and were reviewed by Sreenivasan 
and Stolovitzky [2l| The scale similarity of random fields 
was explored by Novikov [22|, [23|, with a focus on the 
energy dissipation cascade. In Novikov's work, the ratio 
of dissipation averaged over two spheres, one embedded 
within the other, served as a measure of enstrophy par- 
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titioning between larger and smaller scales. The prob- 
ability distribution of these ratios, known as multipliers 
or breakdown coefficients, was shown to relate to multi- 
fractal and statistical measures (moments) of the velocity 
and dissipation fields. A recent review of intermittency 
in multiplicative cascades stresses that this theory is a 
kinematic description and its connection with the real 
dynamics remains unclear [24| . 

Our previous numerical study of particle concentration 
in turbulent flows showed that the particle density field 
is a multifractal on scales comparable to the Kolmogorov 
length scale [l4j. This result suggests that a deeper de- 
scription of the statistical properties of the particle con- 
centration field, based on multiplier PDFs, may also be 
possible. Analytical efforts have suggested that dissipa- 
tion and vorticity in the fluid phase should be locally 
linked with particle concentration [251 ] - Numerical work 
in this regard has demonstrated that preferential concen- 
tration is statistically anticorrelated with low vorticity: 
particles tend to concentrate in regions where enstrophy 
is relatively weak (2|| [2?| • 

In this paper we present a cascade model in the spirit 
of Novikov [22|, [23| that follows the partitioning of pos- 
itive definite scalars associated with both the fluid and 
the particles. Multipliers controlling the partitioning of 
enstrophy and particle density at each step in the cascade 
are drawn from probability distribution functions (PDFs) 
which are determined empirically from direct numerical 
simulations (DNS). Moreover, the multiplier PDFs are 
dependent on, or conditioned by, the particle mass den- 
sity or mass loading. The cascade model then generates 
joint PDFs for particle concentration and enstrophy at 
arbitrary cascade levels. A partitioning correlation prob- 
ability is also applied at each cascade level to account for 
the observed spatial anticorrelation between enstrophy 
and particle concentration [26|, [28| . 

In Section II we describe the cascade model and its 
parameters, which are empirically determined from DNS 
calculations. Details of the DNS equations, and our nu- 
merical methods, are discussed in the Appendix. Results 
are shown in section III, including comparisons of joint 
PDFs of enstrophy and particle concentration as pre- 
dicted by the cascade model with those obtained directly 
from the DNS results. Cascade model PDF predictions 
at Reynolds numbers well beyond the DNS values are 
also presented. In section IV, we summarize our results 
and discuss their implications. 



II. 



CASCADE MODEL 



A turbulent cascade can be envisioned as an hierar- 
chical breakdown of larger eddies into smaller ones that 
halts when the fluid viscosity alone can dissipate eddy ki- 
netic energy. Eddies or similar turbulent structures such 
as vortex tubes are bundles of energy containing vorticity 
and dissipation. These structures start with a size com- 
parable to the integral scale A of the flow, and break down 



in steps to a size comparable to the Kolmogorov scale 
77 before being dissipated away by viscosity. The fluid 
vorticity and dissipation exhibit spatial fluctuations that 
increase in intensity as the spatial scale decreases. This 
phenonemon is known as intermittency and has been ob- 
served in a variety of processes with strong nonlinear in- 
teractions. 

In previous numerical and experimental studies, locally 
averaged intermittent dissipation fields with scale at or 
near r\ were used to quantify the statistical properties 
of multiplier distributions [2l|. Multipliers are random 
variables that govern the partitioning of a positive defi- 
nite scalar as turbulent structures break down along the 
cascade. In these studies the statistical distribution of 
multipliers (their PDF) were shown to be invariant over 
spatial scales that fall within the turbulent inertial range. 
Multifractal properties of the cascading field are deriv- 
able from such multiplier distributions |23fl . and cascade 
models based on the iterative application of multipliers 
to a cascading variable have been shown to mimic inter- 
mittency. 

While invariant with level in the inertial range of a 
cascade, multiplier PDFs might depend on local proper- 
ties of the environment. For instance, Sreenivasan and 
Stolovitzky [2l[ showed that the degree of intermittency 
in dissipation increases with the degree of local strain 
rate, and constructed multiplier distributions for local 
energy dissipation conditioned on the local strain rate. 
The physical mechanism behind this effect is believed to 
be related to vortex stretching dynamics creating intense 
bursts of dissipation. 

All the multiplier PDFs measured by Sreenivasan and 
Stolovitzky 21|, whether conditioned or unconditioned 
by local properties, are well characterized by the (3 dis- 
tribution function, 
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where m is the multiplier variable and /3 is a shape con- 
trolling parameter. A large (3 produces a narrow, delta- 
function-like curve centered at m = 0.5, whereas (3 — 1 
produces a flat distribution between m = and 1. These 
limits for (3 correspond to uniform and highly intermit- 
tent processes respectively. In conditioned multipliers, 
the value of [3 varies with some local property of the 
fluid. 

Concentration of particles in turbulence is a result of 
the active dynamics of eddies on all scales. The process 
depends on the scale of the eddies and the corresponding- 
particle response to those eddies. Intense particle den- 
sity fluctuations, akin to intermittency, were observed 
in a previous numerical study where it was also shown 
that nonuniform particle concentrations have multifrac- 
tal scaling properties [3| . These results strongly suggest 
that a phenomenological cascade model based on multi- 
pliers may adequately describe the particle density field. 
Simulations that have included particle feedback on the 
fluid through the mass loading effect show that damp- 
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ing of local turbulence occurs 0, [2t|. The latter have 
shown that vorticity dynamics is affected locally by par- 
ticle feedback. This interplay between the phases could 
attenuate vortex stretching and, thereby, diminish local 
turbulent intermittency. Multiplier distributions condi- 
tioned on local mass loading should therefore be an inte- 
gral part of a realistic fluid-particle cascade model. 

A. Two-Phase Cascade model 

Below we describe a two-phase cascade model that in- 
corporates simultaneous multiplier processes for parti- 
cle concentration C and fluid enstrophy S, in addition 
to a process that models their spatial anticorrelation. 
The multiplier distributions are conditioned by the local 
particle concentration, as determined empirically from 
DNS fields equilibrated to Re x = 34, 60, 107, and 170. 
The spatial anticorrelation was also quantified from these 
fields. Local measures of particle concentration (C) and 
enstrophy (S 1 ) used are defined in the Appendix. 

A schematic illustration of our two-phase partition- 
ing process is shown in FIG. [TJ The cascading vector 
(5, C) has components representing enstrophy and parti- 
cle concentration. Initially the components are assigned 
the value unity and are associated with a common cell 
having a volume of unity. Each component is partitioned 
into two parts; (msS, (l-mg)S) and {mcC, (l-mc)C), 
respectively, where ms,mc are multipliers for S and C 
whose values are between zero and one inclusive and are 
random members of the corresponding multiplier distri- 
butions. The parts are associated with two daughter cells 
each containing half the volume of the starting cell. In 
the example shown in FIG. [TJ mj and uic are assumed 
to be greater than 0.5. The largest parts of S and C are 
placed in the same daughter cell with probability T (and 
in different cells with probability 1 — T) . This partition- 
ing process is repeated for each daughter cell down the 
cascade until the ratio of the daughter cell size to the 
initial cell size equals a specified cutoff. When this cutoff 
is set to the ratio of the turbulent lengthscales A and 77, 
the cascade corresponds to turbulence characterized by 

i? eA ~ (A/77) 2 / 3 m 
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FIG. 1: Figure depicting the breakdown of a parcel of en- 
strophy (S) and particle concentration (C) into two parcels 
each with half the volume of the parent. The corresponding 
multipliers ras and mc are assumed to be greater than 0.5 in 
this figure. These measures are broken down and distributed 
between the two parcels in one of two ways - the larger por- 
tions are partitioned together with probability F= 0.3 (upper 
figure), or in opposite directions with probability 1 — T= 0.7 
(lower figure). 



Conditioned Multipliers 



The parameters of the cascade model are empirically 
derived from the particle density and enstrophy fields C 
and S as calculated by DNS (see Appendix). The simu- 
lation parameters for four DNS runs representing Re\ = 
36, 60, 104, and 170 are shown in Table [TJ The turbu- 
lence kinetic energy q, the volume averaged dissipation 
e, and A are calculated from the 3-D turbulent energy 
spectrum E{k) and kinematic viscosity v, 



E{k)dk 



(2) 



e = 2v E{k)k 2 dk 
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where k is wavenumber. k max = ^ times the number 
of computational nodes per side is the maximum effective 
wavenumber. Thus k max r\ > 1 indicates an adequate 
resolution of the Kolmogorov scale. 
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Parameter 


Case I 


Case II 


Case III 


Case IV 


0.15 




Nodes/side 


64 


128 


256 


512 






V 


.01 


.003 


.0007 


.0002 


& 01 




Rex 


34 


60. 


104 


170 


0.05 




q 


1.5 


.65 


.28 


.14 







e 

V 


23. 


22.8 


22.4 


23 







1.4 


1.5 


1.45 


1.56 






A 
V 


14.1 


23.3 


45.8 


86.2 


0.15 




r 


.31 


.29 


.27 


.32 


§ 0.1 
y 




D 


.0001 


.00003 


.000007 


.000002 




Vp 


.001 


.0003. 


.00007 


.00002 


^0.05 





TABLE I: Case Parameters for DNS runs. The quantities D 
and v v are defined in the Appendix. Other quantities above 
are defined in Section II. 



The 3-D DNS computational box is uniformly subdi- 
vided into spatial cells 3?7 on a side, and the average value 
of C and S is determined for each cell ( see Appendix ). 
The cells are divided into groups associated with disjoint 
ranges of C. Each cell is then divided into two parts of 
equal volume and averages for C and S are determined 
for each part. The C and S multipliers for each cell are 
evaluated as the ratio of these averages to the averages 
in the parent cell. A conditional multiplier distribution 
p(m) is then determined for each binned value of C from 
the corresponding set of cell multipliers. Plots of p(m) for 
three values of C are shown in FIG.[2j The points repre- 
sent distributions derived from all DNS runs and the solid 
lines are least squares fits to the (3 distribution function 
(Eq. [T|). For the lower values of C, i?eA-independence 
is apparent; only the Re\ = 170 case provided data for 
the largest C range. The plots clearly indicate that the 
intermittency in C is reduced (multiplier PDFs narrow) 
as C is increased. Derived values of (3c(C) and f3s{C) 
are shown as a function of C in FIG. [3l Least squares 
fits to the functional form p 1 exp(p 2 C P3 ) are drawn as 
solid lines and the best fit parameter values for this func- 
tion are tabulated in Table [Til Bounding curves (dashed 
lines) are defined by setting p 2 and p 3 to their 2a lim- 
its, to establish a plausible range of uncertainty in the 
predictions. 



Scalar 


pi 


P2 


pz 


C 


2.7 


.045 


1.02 


S 


9. 


.03 


1.06 



TABLE II: (3 model parameters 

It is certainly of interest that such large solid/gas mass 
loadings as C = 100 appear in the DNS runs at all, given 
published reports that particle mass loading significantly 
dampens turbulent intensity even for mass loadings on 
the order of unity [l|, |J] . These diverse results might be 
reconciled since the particles we study herein are all far 
smaller than the Kolmogorov scale and also have only a 
very small lag velocity relative to the gas. Recall that 
we force the turbulence, as might be the case if it were 
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FIG. 2: Empirically determined conditional multiplier distri- 
butions p(m\C) for particle concentration at three different 
mass loading values, C = 1, 20 and 50. The distributions are 
obtained from bifurcations of cells with a spatial scale equal 
to 377. Results at Rex = 34 ( square ), 60 (triangle), 107 (cir- 
cle) and 170 ( cross ) are overlain. Only the simulation with 
Rex = 170 provided results for C = 50. At each mass loading 
the p(m) at all Reynolds numbers are very well approximated 
with the f3 distribution function ( solid line ). The distribu- 
tion widths narrow as the mass loading increases, indicating 
a decrease in the intermittency. 



being constantly forced by energetic sources operating on 
larger scales than our computational volume. However, 
FIG. [3] strongly suggests an upper limit for C ( ~ 100 ) 
for both (3s and (3c- 

The cascade anticorrelation parameter T was deter- 
mined by counting the number of parent cells within 
which the larger partitions of C and 5* were found to 
share the same daughter cell. This number divided by 
the total number of parent cells defines T. The derived 
r value is approximately constant across the DNS cases, 
as indicated in Table UJ Operationally, the T used in 
the cascade model was determined by taking a simple 
average of the T values in Table [TJ 

Overall, the invariance of V and the (3c (C) and (3s (C) 
functions across our range of Re\ justifies their treat- 
ment as level independent parameters in the two-phase 
cascade model. One caveat remains, which would be of 
interest to address in future work. While it has been 
shown that multiplier distributions leading to (3c and (3s 
are level-invariant over a range of scales within an inertial 
range [21,], our simulations were numerically restricted 
to values of Re in which the inertial range has not yet 
become fully developed. Our reliance on the smallest 
available scales of 3?y to 1.5?y (those providing the largest 
available intermittency) might lead to some concern that 



■5 



1000 - 



100 - 

p 



10 - 




1 

0.001 



0.01 



0.1 







100 



c 



FIG. 3: The (3 parameters as functions of local mass loading 
C for enstrophy and particle concentration at 2>r). Results for 
all DNS cases are indicated as described in FIG. [2] A least 
squares fit of an exponential function to the points over the 
entire mass loading range is shown ( solid line ). Dashed lines 
correspond to the upper and lower limits of the function, and 
are derived using the 2a errors of p2 and P3. 



they were already sampling the dissipation range of our 
calculations, and thus may not be appropriate for a cas- 
cade code. We tested this possibility by calculating mul- 
tipliers for the next largest level bifurcation (677 to 3?y) 
for the Rex = 170 case. The j3 values for those multi- 
plier distributions are slightly larger in value, but consis- 
tent with the C-dependence shown in FIG. [2] (677 scales 
don't provide good distribution functions beyond C ~ 
15). Thus we believe that for the purpose of demonstrat- 
ing this technique, and for the purpose of estimating the 
occurrence statistics of C under particle mass loading, 
our results are satisfactory. For applications requiring 
quantitatively detailed and/or more accurate P{S, C), it 
would certainly be of interest to extend the DNS calcu- 
lations to larger Re, at which a true inertial range might 
be found. 



is P(S, C)ASAC. For quantities varying over orders 
of magnitude, it is convenient to adopt AS — S and 
AC = C, and we will present the results in the form 
P(S, C)SC. 

We started by binning results at spatial scale 3ry, ob- 
tained from the semi-final level of a cascade model run, 
into a uniform logarithmic grid of S, C bins each having 
width A(logS) = A{logC) = 5, with corresponding val- 
ues of AS and AC. The number of 377 cells accumulated 
in each bin was normalized by the total number of such 
cells in the sample to convert it to a fractional volume 
AV(S, C) = P(S, C)ASAC. Then 

AV(S,C) P(S,C)ASAC 



S 2 



A{logS)A(logC) 



P(S,C)SC as <y-»o. 

(5) 

In practice of course, the binning ranges d are not van- 
ishingly small. 

The plots in FIGs. g] and then, show the PDF 
as the volume fraction P(S, C)SC. Cascade levels 9, 12, 
15, and 18 correspond approximately to the Re\ of the 
four simulation cases shown in Table fl] These levels 
were determined from the ratio of A and 77 for each case: 
level = 31og 2 (A/?7). The factor 3 accounts for cascade 
bifurcations of 3D cells, because it takes three partition- 
ings, along three orthogonal planes, to generate eight 
subvolumes of linear dimension one-half that of the par- 
ent volume. That is, 2 is equal to the number of 
r] cells within a 3D volume having linear dimension A 
and (2 levet / 3 ) 2 / 3 is the corresponding Re\. The number 
of cascade realizations is, in turn, equal to the product 
of the number of A-size volumes in the computational 
box and the number of simulation snapshots processed. 
In general it is difficult to generate DNS results with a 
ratio of A and 77 that is an exact power of two. In or- 
der to correctly compare DNS simulations with the cas- 
cade model it was necessary to interpolate between two 
cascade generated P(S, C)SC computed at scale ratios 
(levels) that bracketed the ratios that were actually sim- 
ulated. In FIG. H] we compare iso-probability contours of 
P(S, C)SC predicted by cascade models representing the 
four DNS cases with the same contours derived directly 
from the simulated S and C fields. The agreement is very 
good. 



Predictions at higher Reynolds number 



III. MODEL RESULTS 

The 2D joint probability distribution function or PDF 
of concentration and enstrophy, a fractional volume mea- 
sure, was generated from the cascade model and com- 
pared with results derived directly from numerical DNS 
simulations. The basic probability density P(S, C) gives 
the fractional volume occupied by cells having enstrophy 
S and concentration C, per unit S and C; thus the frac- 
tional volume having C and S in some range AS, AC 



The cascade model was used to generate PDFs at 
deeper levels in order to assess the effect of mass loading 
on the probabilities of high C and S. We generated 256 
realizations of a level 24 cascade, 20 realizations of a level 
30 cascade, and one realization of a level 36 cascade. 

FIG.OJa) shows the average of 256 realizations of a 24 
level cascade, taken to lower probability values. The pro- 
nounced crowding of the contours at the top of the figure 
indicates the effect of particle mass loading on reducing 
the intermittency of C at high values of C. For compar- 
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FIG. 4: Comparisons of cascade model predictions of 
P(S,C)SC with DNS results at 7?e A = 34 (a), 60 (b), 107 
(c) , and 170 (d). Contours indicate probabilities .001, .01, .1 
and .3. Dashed contours are cascade model predictions and 
solid ones are DNS results. 



ison, FIG. (Hb) shows a control run of a 24 level cascade 
with all conditioning turned off. In this control case, the 
exponential tails characterizing intermittent fluctuations 
are seen at both low and high C. 

In order to evaluate the effect of the uncertainties in the 
extrapolations of the j3 curves for C and S on the PDF, 
two cascade runs to level 24 were generated using the pa- 
rameters for the upper and lower dotted curves in FIG. [3] 
In FIG. [6] we show cross-sections of the PDFs produced 
by these runs along the C axis through the distribution 
modes to compare with the same cross-section for a run 
using the nominal parameters in Table [TTJ Both models 
diverge from the mean model beyond C > 40, with the 
upper (lower) curve corresponding to the outside (inside) 
/3c(C) and /3s(C) bounds in FIG. [31 Figure [5] indicates 
that the sensitivity of the PDF to the /3 model parame- 
ters at the 2tr level is only apparent at large C, and all 
models show a sharp dropoff in the probability for C > 
100. 

A crowding effect similar to the one seen in FIG. El^a) 
is shown in FIG. [7] for iso-probability contours equal to 
5 x 10~ 4 , for cascade levels 6, 12, 18, 24, 30 and 36. 

Figures [3Ja) and [SJb) compare ID cuts through the 
modes of the PDFs for cascades of 18 - 36 levels, indi- 
cating that going to deeper levels (higher Re\) results in 
larger intermittency at the low-C end (as expected), re- 
taining the exponential tail characteristic of intermittent 
processes, but the highest particle concentration end of 
the distribution is extended more slowly. Certainly at 
the order of magnitude level, a particle mass loading ra- 
tio of 100 times the gas density appears to be as high as 
preferential concentration can produce. This result could 
be inferred directly from inspection of the conditioned (3 
distributions of FIG. [3) 
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FIG. 5: (a) Cascade model predictions for a 24 level case, 
taken to lower probability levels, using 256 realizations of the 
cascade. Contours are labeled by log(P(S,C)SC). Note the 
crowding of contours at high C values, indicating the high-C 
limit of the process under conditions of mass loading, (b) A 
control cascade to level 24, as in FIG.[5j a ), with conditioning 
turned off. The difference between (a) and (b) clearly shows 
the "choking" effects of particle mass loading on intermittency 
in C. 



IV. SUMMARY 

A two-phase cascade model for enstrophy and parti- 
cle concentration in 3-D, isotropic, fully developed tur- 
bulence with particle loading feedback has been devel- 
oped and tested. Multiplier distributions for enstrophy 
and particle concentration were empirically determined 
from direct numerical simulation fields at Taylor scale 
Reynolds numbers between 34 and 170. These simula- 
tions included 'two-way' coupling between the phases at 
global particle/gas mass loadings equal to unity. The 
shape of all multiplier distributions is well characterized 
by the (3 distribution function, with a value of (3 that 
depends systematically on the local degree of mass load- 
ing. The values of f3 increase monotonically with mass 
loading and begin to rapidly increase at mass loadings 



FIG. 6: ID cuts through the mode of the PDF of FIG. EJa) 
parallel to the C axis, showing the effects of uncertainty in 
the conditioning curve j3c(C). The solid curve is the nominal 
model and the dashed curves are obtained by allowing the 
parameters p2 and pz to take their 2<r extreme values. 
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FIG. 8: (a) ID global cuts through the cascade model PDFs 
P(S, C)SC for runs with 18, 24, 30, and 36 levels, (b) closeup 
of 1-D cuts through high-C regime. 
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FIG. 7: Cascade model predictions for P(S, C)SC = 5 x 1(T 4 
for levels 6, 12, 18, 24, 30, and 36. Contour labels indicate 
the cascade levels. 



greater than 100. 

The C-dependent multiplier distributions were used as 
input to a cascade model that simulates the breakdown, 
or cascade, of enstrophy S and particle concentration C 
from large to small spatial scales. The spatial anticorre- 
lation between enstrophy and particle concentration was 
empirically determined from 3D DNS models and shown 
to be constant with Re\. This constant was used as 
a correlation probability governing the relative spatial 
distribution of S and C at each bifurcation step in the 
cascade model. 

The cascade model we have developed clearly repro- 
duces the statistical distributions and spatial correlations 
observed in our DNS calculations. The cascade parame- 
ter values we have derived appear to be universal within 



the range of Re\ of our simulations. We thus specu- 
late that they can be used to predict approximate joint 
probabilities of enstrophy and particle concentration at 
higher Reynolds numbers, at great savings in computer 
time. For example, a typical DNS run to Re\ = 170 
takes about 170 cpu hours on an Origins 3000 machine, 
while a cascade model to an equivalent level takes 0.1 cpu 
hours. 

We have presented joint probabilites of S and C de- 
rived from cascade runs up to level 36. The contours 
shown in FIG. [Sfa) and FIG. [5] clearly show the effects 
of particle mass loading on the probability distribution 
functions of C in the regimes where C is large. It appears 
that particle mass loadings greater than 100 are rare in 
turbulent flows. 

The properties of the cascade rest on the physics of our 
DNS simulations, and we speculate that two separate ef- 
fects are involved. First, particle mass loading dampens 
fluid motions of all types, decreasing vorticity stretching 
and all other forms of ongoing eddy bifurcation which are 
needed to produce intermittency. Second, as a byproduct 
of this, particle mass loading may alter the Kolmogorov 
timescale locally and shift the most effectively concen- 
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trated particle Stokes number St to a larger value than 
that characterizing particles already lying in the local 
volume, reducing the probability of preferentially con- 
centrating the local particles any further. 
Caveats and Future Work: 

As described in section II, our multiplier distributions 
were taken from the most numerous cells, with the largest 
intermittency, which are at the smallest scales possible 
(furthest from the forcing scale). At Reynolds numbers 
accessible to DNS, a true inertial range is only beginning 
to appear, and while, sampling at the smallest spatial 
scales possible, we are as closely approaching the asymp- 
totic values within the true inertial range as possible, 
where level-independence has been demonstrated in the 
past [21| . it is possible that our values are subject to 
inaccuracy by virtue of being sampled too close to the 
dissipation scale. Any such inaccuracy will affect our 
cascade results quantitatively but not qualitatively. As 
computer power increases, it would be a sensible thing 
to continue experiments like these at higher Re\. 

A more general model that treats enstrophy and strain 
as independent cascading scalars might allow for a higher- 
fidelity particle concentration cascade, since C is known 
to be linked to the difference between these two scalars 
[25l ] (the so-called second invariant tensor II). However, 
such an effort would introduce further complexity of its 
own, as II is no longer positive definite. We consider the 
development of such a model a suitable task for future 
work. 



APPENDIX 

We used an Eulerian scheme developed by Dr. Alan 
Wray to solve the coupled set of fluid/particle equations 
used in this study. This was done to maximize the 
computational efficiency of the calculations and, more 
importantly, to accurately evaluate multipliers over the 
wide range of particle concentrations and enstrophics ex- 
pected. In this study the effects of particle collisions and 
external forces on the particles (e.g., gravity) are not con- 
sidered. The turbulence is spectrally forced at k = \f\A 
such that moments of the Fourier coefficients of the force 
satisfy isotropy up to the fourth order. The instanta- 
neous Navier-Stokes equations describing the conserva- 
tion of mass and momentum for an incompressible fluid 
are 



The compressible equations for the particles are 



V-U = 



<9U 
~dt 



(U- V)U 



VP 

pf 



is\7 2 XJ 



Pf 



(A.l) 



V) (A.2) 



where U is fluid velocity, V is particle velocity, pf and 
p p are the fluid and particle mass densities, v is fluid 
viscosity, P is pressure, and a is the inverse of the particle 
gas drag stopping time r p . 



dpp 
dt 



+ V(p p V)= DV 2 Pp 



(A.3) 



d(p P V) 
dt 



-V(p p VV)= v p V 2 { Pp V)+a Pp {V-V) (A.4) 



where v p is a "particle viscosity" , and D is a "particle 
diffusivity" . The particle diffusivity and viscosity terms 
numerically smooth out particle mass and momentum, 
alleviating the formation of steep gradients of p p that 
can lead to numerical instabilities eg. [311 ]. 

The right hand sides of Eqs. IA.2l and lA.4l contain phase 
coupling terms which are linearly dependent on (U — V). 
The linear form of the coupling follows from the assump- 
tions that the particle size is much less than rj, and that 
the material density of the particles is much greater than 
Pf [2j. Additional contributions to the particle-gas cou- 
plings involving pressure, viscous and Basset forces [29I ] 
have not been added since they are expected to be weak 
in our size regime of interest. The particle field is intro- 
duced with a constant mass density and an initial veloc- 
ity given by the local gas velocity in a field of statisti- 
cally stationary turbulence. All runs are continued until 
the particle statistics (RMS of conentration distribution) 
have equilibrated. 

The particle Stokes number St is defined relative to 
the Kolmogorov time scale r I( as St — t p /t v , and $ = 
Mp/M f is the global mass loading, where M p and M/ are 
the total mass of particles and fluid respectively. In this 
study pf, St, and $ are set to unity, Djv = 0.01, and 
Vpjv = 0.1. Explicitly setting St = 1 guarantees that 
the particles are preferentially concentrated. When $ is 
unity, p p is a surrogate for the local mass loading or local 
concentration factor C . The values of v p and D minimize 
the diluting effects of numerical particle diffusion while 
preventing numerical blowups; their values were deter- 
mined from a set of DNS runs in which their values were 
decreased systematically until numerical instabilities set 
in. 

Eqs. IA.ll - rA~4l are solved using psuedo-spectral meth- 
ods commonly used to solve Naviers- Stokes equations for 
a turbulent fluid. The Fast Fourier Transform (FFT) al- 
gorithm is used to efficiently evaluate the dynamical vari- 
ables U, V and p p on a 3D uniform grid of computional 
nodes with periodic boundary conditions. The computa- 
tional algorithm is parallelized using MPI and is written 
in Fortran 90. All runs for this study were executed on 
SGI Origins supercomputers with up to 1024 processors. 

Enstrophy is defined as 



(A.5) 



where i,j are summed over the three coordinate dimen- 
sions of U. 
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The local spatial average of a scalar over a sample vol- 
ume is estimated as, 

i 

where Fi is the scalar's value on computational node i 
centered within a cube of volume dv and the sum is over 
all n nodes covering the sample volume. We normalized 
this average by the global average value to get a quan- 
tity that measures the scalar's local value relative to its 
mean. In this paper C and S will denote normalized spa- 
tial averages of particle concentration and enstrophy over 
cubes 3?7 on a side. 
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